Get started with fittr

library(fittr)

Why fittr?

In applied research — economics, epidemiology, political science — a typical project involves running dozens of related regressions. You might have several outcome variables, a handful of model specifications (base, extended, with interactions, with additional controls), multiple subsamples, and both weighted and unweighted variants. The combinations add up fast: 3 outcomes x 4 model specs x 2 datasets x 2 weight variants = 48 regressions.

The manual approach is to write out each call individually, or to nest for loops with ad-hoc naming. Both are painful to maintain. When a reviewer asks you to add a covariate to every model, you have to hunt down every call and edit it. It is easy to introduce copy-paste errors, and there is no simple way to preview the full set of runs before committing to them.

fittr takes a declarative approach: you define your model specifications once, list your datasets, and describe each dependent variable in a definition. fittr generates every combination and runs them — optionally in parallel — returning a cleanly named nested list of fitted model objects. Because fittr calls your modelling function with a standard (formula, data, weights) interface, it works with lm(), glm(), mixed-effects models, survival models, GAMs, and anything else that follows that convention.

Building model specifications

The building block in fittr is the model_spec, which represents the right-hand side of a formula. You create one with model_spec():

base <- model_spec(cyl, disp)
base
#> model_spec: ~ cyl + disp

model_spec() accepts unquoted variable names, strings, or R expressions:

model_spec("cyl", "disp")
#> model_spec: ~ cyl + disp
model_spec(I(hp^2), log(wt))
#> model_spec: ~ I(hp^2) + log(wt)

The real power is composing a family of related models from a shared base using pipe-friendly helpers:

models <- list(
  base     = base,
  extended = base |> add_terms(hp, wt),
  reduced  = base |> drop_terms(disp)
)

models
#> $base
#> model_spec: ~ cyl + disp 
#> 
#> $extended
#> model_spec: ~ cyl + disp + hp + wt 
#> 
#> $reduced
#> model_spec: ~ cyl

When you later decide to add a term to base, every derived model inherits the change automatically — no find-and-replace needed.

Other helpers include replace_term() for swapping one term for another, add_interaction() for adding main effects and all their interactions, and update() for formula-style updates:

replace_term(base, "disp", "drat")
#> model_spec: ~ cyl + drat
base |> add_interaction(cyl, disp)
#> model_spec: ~ cyl + disp + cyl:disp
update(base, ~ . + hp - disp)
#> model_spec: ~ cyl + hp

Your first fit

With models defined, you need two more things: a dataset and a list of definitions. Each definition describes one dependent variable — at minimum, the variable name (dv) and the fitting function (fun):

defs <- list(
  mpg = list(dv = "mpg", fun = "lm")
)

results <- fit_models(models, mtcars, defs)

The result is a nested list — outer level keyed by definition name, inner level keyed by model name:

names(results)
#> [1] "mpg"
names(results$mpg)
#> [1] "base"     "extended" "reduced"

Each element is a standard fitted model object, so all familiar R tools work:

coef(results$mpg$base)
#> (Intercept)         cyl        disp 
#> 34.66099474 -1.58727681 -0.02058363
coef(results$mpg$extended)
#> (Intercept)         cyl        disp          hp          wt 
#> 40.82853674 -1.29331972  0.01159924 -0.02053838 -3.85390352

Passing "lm" as fun is a convenient shorthand. fittr wraps it into a function with the signature function(formula, data, weights); a namespaced name such as "fixest::feols" works too. You can also pass a function directly — the only requirement is that it accepts those three arguments.

Multiple dependent variables

To fit models for more than one outcome, add entries to defs:

defs <- list(
  mpg  = list(dv = "mpg",  fun = "lm"),
  qsec = list(dv = "qsec", fun = "lm")
)

results <- fit_models(models, mtcars, defs)
names(results)
#> [1] "mpg"  "qsec"
names(results$qsec)
#> [1] "base"     "extended" "reduced"

Each DV gets the same set of base models. You can extract summaries across models with standard functional tools:

purrr::map_dbl(results$mpg, \(x) summary(x)$r.squared)
#>      base  extended   reduced 
#> 0.7595658 0.8486348 0.7261800

Previewing before fitting

Before committing to a large batch, use preview_fit() to see exactly what will run:

preview_fit(models, mtcars, defs)
#> # A tibble: 6 × 5
#>   def   dv    model    key      formula                    
#>   <chr> <chr> <chr>    <chr>    <chr>                      
#> 1 mpg   mpg   base     base     mpg ~ cyl + disp           
#> 2 mpg   mpg   extended extended mpg ~ cyl + disp + hp + wt 
#> 3 mpg   mpg   reduced  reduced  mpg ~ cyl                  
#> 4 qsec  qsec  base     base     qsec ~ cyl + disp          
#> 5 qsec  qsec  extended extended qsec ~ cyl + disp + hp + wt
#> 6 qsec  qsec  reduced  reduced  qsec ~ cyl

The key column shows the name used to access each result (e.g. results$mpg$base), and formula shows the full formula that will be fitted. This is especially valuable when definitions include extra terms or model modifications — you can verify the formulas are correct before spending time on the actual fits.

Multiple datasets

Pass a named list of data frames to run every model on every dataset:

datasets <- list(
  all   = mtcars,
  light = mtcars[mtcars$wt <= 3, ]
)

results <- fit_models(models, datasets, defs)
names(results$mpg)
#> [1] "base.all"       "base.light"     "extended.all"   "extended.light"
#> [5] "reduced.all"    "reduced.light"

Result keys now include the dataset name: model.dataset. The preview reflects this too:

preview_fit(models, datasets, defs)
#> # A tibble: 12 × 6
#>    def   dv    model    dataset key            formula                    
#>    <chr> <chr> <chr>    <chr>   <chr>          <chr>                      
#>  1 mpg   mpg   base     all     base.all       mpg ~ cyl + disp           
#>  2 mpg   mpg   base     light   base.light     mpg ~ cyl + disp           
#>  3 mpg   mpg   extended all     extended.all   mpg ~ cyl + disp + hp + wt 
#>  4 mpg   mpg   extended light   extended.light mpg ~ cyl + disp + hp + wt 
#>  5 mpg   mpg   reduced  all     reduced.all    mpg ~ cyl                  
#>  6 mpg   mpg   reduced  light   reduced.light  mpg ~ cyl                  
#>  7 qsec  qsec  base     all     base.all       qsec ~ cyl + disp          
#>  8 qsec  qsec  base     light   base.light     qsec ~ cyl + disp          
#>  9 qsec  qsec  extended all     extended.all   qsec ~ cyl + disp + hp + wt
#> 10 qsec  qsec  extended light   extended.light qsec ~ cyl + disp + hp + wt
#> 11 qsec  qsec  reduced  all     reduced.all    qsec ~ cyl                 
#> 12 qsec  qsec  reduced  light   reduced.light  qsec ~ cyl

Weighted regressions

Supply weight_var to run both weighted and unweighted versions of every model:

results <- fit_models(models, mtcars, defs, weight_var = "wt")
names(results$mpg)
#> [1] "base.no_wt"     "base.wt"        "extended.no_wt" "extended.wt"   
#> [5] "reduced.no_wt"  "reduced.wt"

Each model runs twice — once unweighted (.no_wt) and once with the weight variable (.wt). With multiple datasets, the key becomes model.dataset.wt.

If you only want the weighted variants, set weighted_only = TRUE:

results <- fit_models(models, mtcars, defs, weight_var = "wt", weighted_only = TRUE)
names(results$mpg)
#> [1] "base.wt"     "extended.wt" "reduced.wt"

Customising definitions

So far, every DV has used the same set of base models. Definitions support three optional fields that let you customise models per-DV without changing the shared base.

Extra terms

extra_terms adds terms to every model for a specific DV — useful for DV-specific controls:

defs <- list(
  mpg  = list(dv = "mpg",  fun = "lm", extra_terms = c("am", "gear")),
  qsec = list(dv = "qsec", fun = "lm")
)

preview_fit(models, mtcars, defs)
#> # A tibble: 6 × 5
#>   def   dv    model    key      formula                               
#>   <chr> <chr> <chr>    <chr>    <chr>                                 
#> 1 mpg   mpg   base     base     mpg ~ cyl + disp + am + gear          
#> 2 mpg   mpg   extended extended mpg ~ cyl + disp + hp + wt + am + gear
#> 3 mpg   mpg   reduced  reduced  mpg ~ cyl + am + gear                 
#> 4 qsec  qsec  base     base     qsec ~ cyl + disp                     
#> 5 qsec  qsec  extended extended qsec ~ cyl + disp + hp + wt           
#> 6 qsec  qsec  reduced  reduced  qsec ~ cyl

Notice that am and gear appear in every formula for mpg, but qsec’s formulas are unchanged.

Additional models

additional_models adds models that only apply to a specific DV:

defs <- list(
  mpg  = list(
    dv = "mpg", fun = "lm",
    additional_models = list(
      kitchen_sink = model_spec(cyl, disp, hp, wt, drat, am, gear)
    )
  ),
  qsec = list(dv = "qsec", fun = "lm")
)

preview_fit(models, mtcars, defs)
#> # A tibble: 7 × 5
#>   def   dv    model        key          formula                                 
#>   <chr> <chr> <chr>        <chr>        <chr>                                   
#> 1 mpg   mpg   base         base         mpg ~ cyl + disp                        
#> 2 mpg   mpg   extended     extended     mpg ~ cyl + disp + hp + wt              
#> 3 mpg   mpg   kitchen_sink kitchen_sink mpg ~ cyl + disp + hp + wt + drat + am …
#> 4 mpg   mpg   reduced      reduced      mpg ~ cyl                               
#> 5 qsec  qsec  base         base         qsec ~ cyl + disp                       
#> 6 qsec  qsec  extended     extended     qsec ~ cyl + disp + hp + wt             
#> 7 qsec  qsec  reduced      reduced      qsec ~ cyl

The kitchen_sink model appears only under mpg.

Model modifications

model_mods lets you modify specific base models for one DV only. Each entry is a function that takes a model_spec and returns a modified one:

defs <- list(
  mpg  = list(
    dv = "mpg", fun = "lm",
    model_mods = list(
      extended = \(m) drop_terms(m, hp)
    )
  ),
  qsec = list(dv = "qsec", fun = "lm")
)

preview_fit(models, mtcars, defs)
#> # A tibble: 6 × 5
#>   def   dv    model    key      formula                    
#>   <chr> <chr> <chr>    <chr>    <chr>                      
#> 1 mpg   mpg   base     base     mpg ~ cyl + disp           
#> 2 mpg   mpg   extended extended mpg ~ cyl + disp + wt      
#> 3 mpg   mpg   reduced  reduced  mpg ~ cyl                  
#> 4 qsec  qsec  base     base     qsec ~ cyl + disp          
#> 5 qsec  qsec  extended extended qsec ~ cyl + disp + hp + wt
#> 6 qsec  qsec  reduced  reduced  qsec ~ cyl

The extended model for mpg no longer includes hp, but the extended model for qsec is unaffected.

Model-specific regression functions

When different models need different estimators, pass fun as a named list. Names correspond to model names, and .default serves as the fallback:

models <- list(
  ols   = model_spec(cyl, disp),
  logit = model_spec(cyl, disp)
)

defs <- list(
  am = list(
    dv  = "am",
    fun = list(
      logit    = \(f, d, w) glm(f, d, family = binomial, weights = w),
      .default = \(f, d, w) lm(f, d, weights = w)
    )
  )
)

results <- fit_models(models, mtcars, defs)
class(results$am$ols)
#> [1] "lm"
class(results$am$logit)
#> [1] "glm" "lm"

Models that need more than a formula

Some models are not fully described by the right-hand side of their formula. ordinal::clm() takes a second formula in its nominal argument; fixest::feols() puts fixed effects after a |. A model_spec can carry both, so these models can be composed, modified per DV and previewed like any other.

Extra arguments

set_args() attaches arguments that fittr passes to the fitting function along with the formula, data and weights. With it, the model-specific functions example above doesn’t need closures:

models <- list(
  ols   = model_spec(cyl, disp),
  logit = model_spec(cyl, disp) |> set_args(family = "binomial")
)

defs <- list(am = list(dv = "am", fun = "glm"))

preview_fit(models, mtcars, defs)
#> # A tibble: 2 × 6
#>   def   dv    model key   formula         args                   
#>   <chr> <chr> <chr> <chr> <chr>           <chr>                  
#> 1 am    am    logit logit am ~ cyl + disp "family = \"binomial\""
#> 2 am    am    ols   ols   am ~ cyl + disp  NA

The args column shows each model’s arguments. They survive add_terms() and the other helpers, and can be changed per DV with model_mods, using set_args(), drop_args() or modify_arg().

An argument given as a model_spec is passed as a one-sided formula, so a formula-valued argument can be built from terms and changed with the same helpers. For example, an ordinal model with nominal effects:

base <- model_spec(age, income, educ)

models <- list(
  base    = base,
  nominal = base |> set_args(nominal = model_spec(educ))
)

defs <- list(
  trust = list(
    dv = "trust", fun = "ordinal::clm",
    # For this DV, age also has a nominal effect
    model_mods = list(nominal = \(m) modify_arg(m, "nominal", \(a) add_terms(a, age)))
  )
)

A fitting function you write yourself receives the arguments too, so it needs ... (or arguments of the same names): function(f, d, w, ...).

Arguments that a modelling function evaluates inside the data, notably subset and offset, do not work this way. Use offset(x) as a term instead, and pass a separate dataset in place of subset.

Formula parts

Many packages split the right-hand side into sections with |: fixest’s fixed effects and instruments, betareg’s precision model, pscl’s zero-inflation model. Add these with add_part() or set_parts(). Parts are placed after the main terms, in order; what each one means is up to the modelling function.

library(fixest)

base <- model_spec(hp, wt)

models <- list(
  pooled = base,
  fe     = base |> add_part(model_spec(cyl)),
  fe2    = base |> add_part(model_spec(cyl, gear))
)
# Cluster standard errors by cylinder in every model
models <- lapply(models, set_args, cluster = model_spec(cyl))

defs <- list(mpg = list(dv = "mpg", fun = "feols"))

preview_fit(models, mtcars, defs)
#> # A tibble: 3 × 6
#>   def   dv    model  key    formula                    args           
#>   <chr> <chr> <chr>  <chr>  <chr>                      <chr>          
#> 1 mpg   mpg   fe     fe     mpg ~ hp + wt | cyl        cluster = ~ cyl
#> 2 mpg   mpg   fe2    fe2    mpg ~ hp + wt | cyl + gear cluster = ~ cyl
#> 3 mpg   mpg   pooled pooled mpg ~ hp + wt              cluster = ~ cyl

results <- fit_models(models, mtcars, defs)
etable(results$mpg)
#>                               fe              fe2            pooled
#> Dependent Var.:              mpg              mpg               mpg
#>                                                                    
#> hp              -0.0231 (0.0138) -0.0342 (0.0217) -0.0318* (0.0052)
#> wt                -3.181 (1.167) -2.792. (0.9309)  -3.878* (0.6999)
#> Constant                                            37.23** (3.061)
#> Fixed-Effects:  ---------------- ---------------- -----------------
#> cyl                          Yes              Yes                No
#> gear                          No              Yes                No
#> _______________ ________________ ________________ _________________
#> S.E.: Clustered          by: cyl          by: cyl           by: cyl
#> Observations                  32               32                32
#> R2                       0.85722          0.86376           0.82679
#> Within R2                0.46632          0.47643                --
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The term helpers change the main terms and leave the parts alone. (On a plain formula, update() would wrap the whole right-hand side in parentheses, merging the sections.) To change a part given as a model_spec, apply the same helpers to it through modify_part():

model_spec(hp, wt) |>
  add_part(model_spec(cyl, gear)) |>
  add_terms(drat) |>
  modify_part(1, \(p) drop_terms(p, gear))
#> model_spec: ~ hp + wt + drat | cyl

A part can also be a formula or a string, e.g. fixest’s instrumental variables section, add_part(wt ~ drat).

You still need a closure for anything that depends on the data itself, such as a survey design built from d:

fun = \(f, d, w, ...) {
  design <- survey::svydesign(ids = ~ 1, weights = w, data = d)
  survey::svyglm(f, design = design, ...)
}

Parallel execution

When you have many models to fit, use the workers argument to parallelise:

results <- fit_models(models, mtcars, defs, workers = 4)

This sets up a future::multisession plan internally and restores the previous plan on exit. You can choose a different strategy via the plan argument:

results <- fit_models(models, mtcars, defs, workers = 4, plan = future::multicore)

Alternatively, set up your own plan before calling fit_models() and leave workers at its default (NULL). In this mode, each DV is processed sequentially while models within each DV run in parallel:

future::plan(future::multisession, workers = 4)
results <- fit_models(models, mtcars, defs)
future::plan(future::sequential)

Working with results

Because fittr returns standard model objects, all familiar R tools work on the results:

# Redefine models and defs for this section
models <- list(
  base     = model_spec(cyl, disp),
  extended = model_spec(cyl, disp, hp, wt)
)

defs <- list(
  mpg = list(dv = "mpg", fun = "lm")
)

results <- fit_models(models, mtcars, defs)

Extract coefficients across all models for a DV:

purrr::map(results$mpg, coef)
#> $base
#> (Intercept)         cyl        disp 
#> 34.66099474 -1.58727681 -0.02058363 
#> 
#> $extended
#> (Intercept)         cyl        disp          hp          wt 
#> 40.82853674 -1.29331972  0.01159924 -0.02053838 -3.85390352

Compare R-squared values:

purrr::map_dbl(results$mpg, \(x) summary(x)$r.squared)
#>      base  extended 
#> 0.7595658 0.8486348

Or use any package that works with model objects — confint(), predict(), broom::tidy(), modelsummary, and so on.

Next steps

This vignette covered the core workflow. fittr also supports:

• Complex dependent variables like "log(mpg)" or "Surv(time, event)" — pass them as strings in the dv field.
• Complex predictor terms like s(x) (GAM smooths), (1 | group) (random effects), or I(x^2) — pass them directly to model_spec().
• Suppressing the intercept with model_spec(..., intercept = FALSE).
• Programmatic term injection using !! inside model_spec().

See the function documentation for fit_models(), model_spec(), and preview_fit() for full details.